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Abstract. We present a theoretical study of the equilibrium ordering in a 3D XY 
nematic system with quenched random disorder. Within this model, treated with the 
Replica trick and Gaussian variational method, the correlation length is obtained as 
a function of local nematic order parameter Q and the effective disorder strength T. 
These results, £ ~ Q 2 e 1 ^ and £ ~ (l/r)e~ r , clarify what happens in the limiting 
cases of diminishing Q and T, that is near a phase transition of a pure system. In 
particular, it is found that quenched disorder is irrelevant as Q — > and hence does 
not change the character of the continuous XY nematic-isotropic phase transition. We 
discuss how these results compare with experiments and simulations. 
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1. Introduction 

Classical nematics possess long range orientational order as a result of anisotropically 
biased pair interaction between their long rod-like molecules. Below a certain 
temperature such systems are in the ordered phase, where the rods are on average aligned 
with their long axes parallel to each other. Macroscopically this preferred direction is 
defined by a unit vector called the director, labelled as n [T]. In this paper we investigate 
the effects a quenched random field - disorder coupled to the local order parameter - 
has on the ordering of XY nematics. This is a special sub-class of systems where the 
orientational ordering direction is confined to a plane (hence the name XY), while the 
overall system dimensionality may remain the usual 3D (or any other). This is of course 
fully analogous to the famous XY magnets E] , but with the quadrupolar (nematic) 
order rather than a dipolar. Physically such systems are increasingly common in thin 
planar cells of a liquid crystal, or in new liquid crystalline elastomers made into thin 
free standing films with the director confined in their plane. 

The presence of quenched disorder is implicit in many systems, as different as 
random ferromagnets |T]. vortex latices in type-// superconductors [3] and randomly 
crosslinked nematic elastomers 6J. It arises from impurities whose distribution in the 
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sample and the constraints they impose on the local field conformation do not change 
on the time scale relevant for the evolution of the main ordering. 

The influence of quenched disorder is well understood from a theoretical point of 
view. In the limit of weak disorder! , it was first shown by Larkin that impurities lead to 
the destruction of long range order [Sj. This result was then generalized by Imry and Ma 
who demonstrated that, in less than four dimensions, an arbitrarily weak random field 
confines the long range order to hold only at length scales smaller than the "domain size" 
£. The system breaks into uncorrelated regions and therefore loses overall order and 
alignment j3J[7]. Much work has been done on XY random anisotropy magnets, showing 
magnetization correlations to decay exponentially [HE]: ( m (0) • m(r)) ~ e~ r ^. More 
sophisticated analysis was done on flux lattices with a two component order parameter 
(N = 2, d = 3) whose Hamiltonian incorporates an elastic term and a random potential 

It shows that correlations in fact decay as a power- law giving rise to quasi- 

long range order (QLRO). A functional renormalization group analysis gives a similar 
behavior for XY ferromagnets, but also finds that quasi-long range order is absent in the 
case of a magnetization field with more than three components ^U]- Our work follows 
a similar line to Ref [H] and also predicts QLRO. However, another non-pertubative 
functional renormalization group analysis disputes the presence of QLRO in iV = 2 
d = 3, where the system is found to always stay disordered |TT]. The authors find that 
for N = 2 the critical dimensionality below which QLRO exists to be 3.8. 

In many systems with such glassy nature of macroscopic ordering, the experimental 
observation of key theoretical predictions is difficult and often indirect. The loss of 
the long range orientational order is directly seen in equilibrium polydomain nematic 
elastomers. These materials, formed by crosslinking liquid crystalline polymers in their 
isotropic phase, have inherent quenched disorder. There are two plausible explanations 
why disorder is present, related to each other. Residual heterogeneous strains established 
at the moment of crosslinking give rise to random stresses (a), acting locally on the 
nematic order [T2~] . Alternatively (b), the locally anisotropic crosslinks provide randomly 
oriented anisotropy axes whose disordering effect competes with Frank elasticity, which 
favors a uniform director configuration j^J. There is also a model essentially treating 
an ad hoc random temperature in a nematic system ^Hj but its approach, neglecting 
the fundamental couplings of nematic ordering to the underlying elastic network is not 
relevant for our work. It may be possible that a careful analysis of director correlations 
could resolve the physical origin of quenched disorder due to the network crosslinks, 
since computer simulations based on model (a) predict a faster than exponential decay 
of director correlations [12] , whereas simulations based on model (b) find the decay to be 
exponential ^3]. Nevertheless, in both simulations the disorder correlation length has 
been a rapidly decreasing function of the effective disorder strength. The dependance is 
exponential and therefore poses the problem that the correlation length does not diverge 

! The empirical criterion of "weak disorder" is when the quenched impurities can be treated as acting 
on the background of establishing mean ordering field - as opposed to the strong disorder when one 
cannot define the usual order parameter in the same way as in a pure system. 
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at the limit of vanishing disorder. 

In practice, the polydomain nematic elastomers have the following features. On 
cooling below the nematic-to-isotropic transition temperature (Tjvj) a high order 
parameter phase is formed. Measurement of the order parameter is difficult, with 
perhaps only NMR a suitable technique fSl QHj- There are differing reports in the 
literature over the years, in all cases finding the transition to be continuous - approaching 
a critical point in better quality samples. It is well established that deep in the nematic 
phase the orientation of the nematic director is very non-uniform: n varies continuously 
across the sample, pointing roughly along one direction across a very small region of 
space: following a classical spin-glass pattern which is perhaps misleadingly called a 
"polydomain" structure. Correlations between director orientations at different points 
in space decay rapidly and eventually vanish at distances much larger than £ so that 
long range order is eventually lost This characteristic length scale £, is often called 
the domain size or correlation length. It is well established to be of the order of microns 
|18| ITTJ] in many different systems. Therefore, light passing through such a sample 
is strongly scattered on birefringent regions with randomly oriented optical axis. As 
a result, such a sample is completely transparent at high temperatures, but becomes 
opaque below T NI due to the multiple scattering of the disordered nematic phase. The 
dependence of £ on the magnitude of the order parameter (Q) and the effective strength 
of quenched disorder (T) is the main focus this work. The experimental work quoted 
above was done in thin films leading us to consider directors that are confined to the 
x?/-plane, though being dependent an all three spatial coordinates. This choice would 
also adequately describe nematic membranes [20J |2T] • 

The structure of the paper is as follows. In Section 2 we summarize a physical model 
of quenched disorder in nematic systems following |17j . Section 3 applies the Gaussian 
Variational Method [22J to this problem. We find that the Replica symmetry is broken in 
this system and most of the calculations are performed within the Hierarchical Replica 
Symmetry Breaking framework. In Section 4 we solve the variational equations and 
then obtain £ as a function of Q and T. We also examine how the phase transition 
is changed by presenting a predicted Q(T) plot. Finally, in Section 6, we conclude by 
discussing our results and comparing them with experiments. 

2. Model 

2.1. Sources of quenched disorder 

In nematic elastomers crosslinked in the isotropic phase the network crosslinks act as 
sources of quenched orientational disorder. The crosslinks always contain anisotropic 
groups that provide easy anisotropy axes k, similar to the sources in random anisotropy 
magnets We assume it is favorable for the local director to align along k in the 
vicinity of the crosslink. The orientation of the anisotropy well as the spacial 

distribution of the crosslinks are quenched variables. They remain unaltered with time 
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Figure 1. Schematic representation of how crosslinks provide easy anisotropy axes 
{k}. The nematic director is forced to be aligned, in the vicinity of the crosslink, with 
the axes, which are represent by the arrows. Both the orientation of {k} as well as the 
positions of the crosslinks {Ri} are random. Since the crosslinks are confined by the 
network topology, they add quenched disorder to the nematic system. 



even if external conditions, such as temperature, change. 

We follow previous work ^7J modeling the local coupling between the nematic 
order and the random field. For a crosslink positioned at Rj with anisotropy axis kj a 
free energy = —7 ki • Q ■ kj is attributed to it, where 7 is the coupling strength and 
Qij = Q (niUj — I 5ij) is the XY nematic tensor order parameter at spatial position 
Ri. The magnitude Q of the nematic order parameter is different from the Edwards- 
Anderson order parameter of spin glasses, but serves the purpose: Q — 1 refers to a 
completely ordered state whereas Q = to a disordered one [?]. Employing a continuum 
density of impurities p (r) = J2r, — ^0 an d substituting the full expression of 
we get the coarse-grained form of random field energy for the whole sample: 

Fr -f- = -7^k i -g-k i = - J d 3 r 7 Qp(r)(k-n) 2 + l -N cvl Q, (1) 

where N CT is the total number of impurities (crosslinks). The positive linear term 
(iV cr 7 Q) is a byproduct of the requirement that Qij is traceless. It cancels out when 
random disorder is treated properly in Subsection \2.2\ and it is dropped for now. 

When the (one-constant approximation) Frank elasticity term is included, the full 
continuum Hamiltonian of the nematic system reads: 



H = d s r 



l -K{Vnf - 7 Qp(r)(k-n) 2 



(2) 



A simple dimensional argument gives K ~ ksT/a, where a is the interatomic 
distance below which the continuum limit of elasticity no longer holds (essentially, the 
nematic coherence length pQ). Both the microscopic and the phenomenological Landau- 
DeGennes theories of nematic transition give the elastic constant K to scale as Q 2 for 
Q -C 1. Combining these two estimates, we take that for small Q the elastic constant 
is approximated by K ~ k^TQ 2 /a. Finally, it is noted that the local magnitude of 
the order parameter is taken to be constant across the sample and VQ = 0, that is, 
the effect of quenched disorder lies in the resulting equilibrium director texture. It is 
assumed that the sample is chemically homogeneous and boundary effects are neglected. 
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This Hamiltonian has been shown JB] to explain domain formation in nematic 
elastomers using the classic Imry-Ma argument [7j. Even weak random field destroys 
long range order in length scales greater than the correlation length ~ K 2 / pq(^Q) 2 , 
where in our notation p is the number of crosslinks per unit volume. It is a well known 
fact that the correlation length of the order parameter fluctuations should diverge close 
to the phase transition. In direct contrast, when the Imry-Ma argument is applied to 
nematic elastomers the disorder correlation length tends to zero as the order parameter 
vanishes (£<j oc Q 2 ) because K is a quadratic functions of Q. In Section HI we are able 
to clarify this rare shortcoming of the Imry-Ma argument. 



2.2. Replica method 



We are interested in results averaged over the random distributions of the quenched 
variables p(r) and k because these variables cannot be controlled experimentally. In 
other words, one looks for the macroscopic averages over a variety of realizations of 
p(r) and k. Crosslinking the sample above T^i makes the easy anisotropy axes point 
at random directions in the XY plane, with an isotropic probability of orientation 
P(k) = j-. The crosslinks are dispersed randomly in the sample with the density p(r) 
and the probability that a particular distribution p occurs Gaussian [ 

(Ap) 2 



P[p] 



exp 



- / d 6 r- 



2p 



(3) 



where Ap = p — p is the deviation of local density from p , the mean density of 
crosslinks. 

Averaging the logarithm of the partition function Z to obtain the free energy is not 
algebraically possible, so the Replica Trick [22] is employed to overcome this difficulty 
The expression for the free energy arising from disorder then reads: 



k B T (log Z) p>k 



-k B T^- (Z m ) 
Om 



m=0 



j rp 9 

om 



m=0 



m „ 

] [ / Dn a exp [-/3H rep \ , 



(4) 



where we now have m identical "replicas" of the system, labelled by the index a. 
Equation (0J) provides the definition of the Replica Hamiltonian, which no longer depends 
on k and p. After averaging over disorder this Hamiltonian is found to be: 

1 



rep 



= E 



,6=1 



d d r 



-K (Vn a ) 2 5 ab - T (n a • n b f 



(5) 



where subscribes a and b are the replica indexes and m is the number of replicas that 
is set to zero at the end of the calculation. Parameter T, arising from completing 
the Gaussian square with the random-field term in eq. (j2J), reflects the strength of the 
disorder and has a quadratic dependance on the order parameter magnitude: 



8k B T 



(6) 
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It is noted that all replicas are assumed to have equal disorder strength and equal 
magnitude of the order parameter; i.e. 7 a = % and Q a = Qb for all a and b. 

The simplest way to ensure that the director is a unit vector confined in a plane 
is to parameterize it by a single angle 9: n = {cos6(r) , sin#(r)}. In this notation the 
Replica Hamiltonian, now a functional of 9, reads: 



m „ 

H rep [9(r)] = J2jd 3 r 

a,b 



K o 



ab 



r 

— cos2( 
2 v 



(7) 



where the irrelevant constant term — TV/2 arising from the conversion of cos 2 9 to cos 29 
is dropped. 



3. Gaussian Variational Method 

The cosine in the random field term of Eq. (J2J) posses an obstacle to the development of 
the model. The correlation function and the free energy can be obtained from Eq. (@J) 
as long as the Replica Hamiltonian is quadratic in 9. This is clearly not the case in 
our model and to overcome this obstacle we employ the Gaussian Variational Method 
(GVM), first used in the context of random manifolds by Mezard and Parisi [2*2*] . 



3.1. GVM in Nematic Systems 



A model Hamiltonian Hq is assumed to describe the system, meaning that all its physical 
properties, such as correlations and free energy, are fully defined by H [9] instead of the 
original H rep [9]. It is essential that Hq is quadratic in the fields 9(q), in Fourier space: 

d 3 q 

2 — / (2-j : 



(3H 



a.b 



9 a (q)G-J(q)9b{q). 



(8) 



Here and throughout the paper we use the shorthand notation (3 = l/ksT that makes 
all energy functions non-dimensional. The (unknown) propagator can be written as: 

G-J = f3Kq 2 5 ab - a ab , (9) 

where K q 2 5 a b is chosen to match the elastic contribution in H rep , and a a b is a set of yet 
undetermined parameters approximating the random field effects. The true free energy 
of the system can be cast as an expansion in {H rep — Hq), assumed small: 

F w F + (H rep -H ) Ho , (10) 

where F = — ksT (det G) 1 ^ 2 and the angular brackets indicate averaging with the 
Gaussian e _/3H ° as measure. The ultimate aim is to minimize F[G] with respect to 
variations in G a b and obtain the best upper bound for F. Using Eq. (jSJ), the average of 
Eq. (|TU|) leaves the variational energy density: 



f3F v , 



trlogG + ^/?Kg 2 G 



I m r r 

-J2p t ex p - 2 / - 2G - b + Gbh ) 
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where notation J is a shorthand for the full J d 3 q/(2ir) 3 and the irrelevant constant 
(Hq) = —JvbT/2 has been dropped. Before proceeding with the minimization it is 
convenient to define a new function that appears regularly in the random field term: 

B a b = / (G aa — 2G a .b + Gbb) ■ (11) 

Jq 

It is related to the correlations between replicas since ([9 a (r) — b {r)] 2 ) = B ab . 

Demanding F var to be stationary with respect to variations on G a b provides a 
relation that determines a a b- 

5PF VI 



var 



5G 



cd 



G- b 1 = f3Kq 2 5ab-W l e - 2B ^-5abY, e ~ 2Bac ) ( 12 ) 



Comparing this expression with F var , it is easily seen that the stationary equation for 
the matrix of parameters a takes the form: 



(Tab = W 



e" 2B - - 5 ab £ 

c 



(13) 



The above equation has a reassuring property, that a is zero when there is no disorder 
in the system (r = 0). After all, a a b was introduced to approximate the random field 
imposed by the crosslinks. 

It is convenient to decompose the equation into two separate conditions 22]; one for 
the off-diagonal and the separate for the purely diagonal part of a a b- The off-diagonal 
part is easily obtained by setting the Kronecker delta to zero in Eq. (fT3*j) : 

o-a^b = 4{3TVe- 2B « b . (14) 

This is a transcendental equation because B a b is itself a function of G a b and therefore of 
a a b as well. To be able solve it a certain form of a a b must be assumed so that G~^ can 
be inverted to obtain B ao . The diagonal part is fixed by noting that summation over the 
free index a in Eq. (JT3j) is zero because B a b = Bb a . Hence after finding the off-diagonal 

elements the condition 

m 

J2°ab = 0. (15) 
a=l 

provides a way to determine the diagonal matrix element a aa - 



3.2. Replica symmetry 

In a first attempt to determine the variational parameter a a b from the stationary 
equation, let us assume all the non-diagonal elements to be equal to a a ^b = cr, a constant 
which will be determined by Eq. (I14j) . Using Eq. ()15|) we then find the diagonal part to 
be: a aa = (1 — m) a . This scheme based on the constant-a assumption is called the 
Replica symmetry (RS) limit. Under this assumption, the model Hamiltonian acquires 
the simple form: 

Gab(q) = (Kq 2 + ma) S ab (16) 
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where K = (3K and l ao is the matrix with all its elements equal to one. The matrix 
can be trivially inverted: 

Gab = ~Kf 5ab + Wq~ A lab ^ 
The stationary equation (fT^ is easily solved because B ao does not depend on a. This is a 
unique feature of RS in this problem. We shall see later that less symmetric forms of a ab 
generate a transcendental equation more difficult to deal with. The replica symmetric 
solution is: 

a-Wev^l^-VT (18) 

where q ma x — 2n/a is the ultraviolet cut-off in momentum space. If we now use the 
earlier estimate for the Frank elastic constant, K ~ ksTQ 2 /a, the solution becomes 

a = 4/3TV exp f-^J °c Q 2 exp[-4/vrQ 2 ] , (19) 

where the quadratic order parameter dependance of T has been also substituted. This 
function has a singular dependance on the order parameter, dropping to zero as Q — > 0. 
As a result, the disorder is irrelevant close to the nematic-isotropic transition. 

A necessary condition for the RS solution to be applicable is that the disorder energy 
must be stable for infinitesimal variations from that solution. The stability analysis was 
first introduced in the context of spin glasses [2H] and later applied to the GVM [22J. 
First, the variational parameter a ab is allowed to vary: a ao = a^jf + e a b, with the replica 
symmetric ansatz labelled as o~ RS and e a b the infinitesimal deviation from it. F var (o~) is 
then expanded to second order in e ao , naturally, without the linear term, represented by 
the zero in the Eq. ()12j) : 

F var = F var (e = 0) + - / ^2 H abtdc (q) e ab e cd . 

^1 a,b,c,d 

The term second-order in e ab involves the four-dimensional tensor of coefficients H ao c i c , 
which is called Hessian and plays a vital role in determining the stability of the this 
solution: for RS to be stable, this term must be positive definite for arbitrary e. This 
is ensured as long as all the eigenvalues of the Hessian are non-negative. In our XY 
system there are only two non-trivial eigenvalues. One of them, the so-called replicon 
eigenvalue, diverges to negative infinite [2~2] |9~]: 

\ _ 1 _ f rf3 g 1 

° J (2tt)2 (Rq*y 

Therefore the replica symmetric solution does not give an established free energy 
minimum and is not appropriate. 

3.3. Hierarchical RSB 



The stability analysis of the RS solution, as well as much previous work on disorder 
systems, invites us to look for the hierarchical replica symmetry breaking solutions [27]. 
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This is known to be the real equilibrium state (i.e. has both Hessian eigenvalues positive 
definite) of spin glasses [2B] and random manifolds It is reasonable to expect 

a similar behavior in our model and therefore the hierarchical method was chosen to 
explore replica symmetry breaking. 

In this model the matrix a a b is assumed to have a nested block-diagonal form with 
off- and on-diagonal parts map to cr(v) and a respectively, with a continuous index 
variable v G [0,1]. As a result and its inverse also have this hierarchical form 
and the quantity of interest when solving the stationary Eq. (fTf|) is now written as 
B ab = 2 J q [g-g(v)\, where 

2 



~9~9(v) = w - 9 — + / dua\u) — , (20) 

K q 2 - a + cr(l) J„ V-^<r - (J + (a) + [a\(v) J 

where the shorthand notations a' = da/dv and [o~](v) = va{y) — dua{u) are used, 
after Parisi Finally the condition on the diagonal part of o in Eq. ()15|) becomes 

dva(v) = (a). (21) 



o 



It should always be kept in mind that the variable v determines which diagonal block 
an element a a b belongs to. The smaller it is, the closer the element is positioned to 
the diagonal. More interestingly, we shall see later that the long range behavior of the 
system is specifically associated with these small values of v. 

The original stationary equation is written in terms of hierarchical matrices as: 

a{v) = WV exp {-4 jT[fif(?) - <^)]} • (22) 

Substituting Eqs. (gDj) in 



a{v) = 4pTV exp |~ 4 ^ 



1 + f 1 °M du 1 \ (23 ) 

Kg* - a + a(l)^ J v [Kg* - a + <<7> + [*}(u)} 2 \ J ' 1 ! 

Parameter v appears only on the lower limit of the w-integral in Eq. (|23|) . As a result 
the stationary equation is greatly simplified by differentiating it with respect to v; 

= 4cT ^ / \k 2~Tfl7 ,12 . ( 2 4) 
J q [A q 2 + [a]{v)} 2 

where it is noted that Eq. (|2*2*|) shows that a(v ) has no g-dependance. This stationary 
equation for a{y ) has two solutions. The first one is obvious: o~'{y) = 0. Had this 
been the real relevant solution, the unique RS would have been the only possibility (the 
parameter a does not depend on v). However, a second solution also exists for o'{v ) ^ 0, 
and it is given by: 

1 = 4<t(t;) [ ^ q - . => a(v) = 2nK 3 / 2 [a]^ 2 (25) 

V } J 2n 2 [Rq 2 + [a](v)} 2 K ' M 1 ' 

Differentiating Eq. (|2*H)l w.r.t. v and observing that [cr]'(f) = va'(v), gives the second 
solution for a: 

[a]{v) = n 2 K 3 v 2 => a(v) = 2n 2 K 3 v . (26) 
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Figure 2. The function a(v), plotted for two different disorder strengths (dashed line 
for higher T). The weakening disorder decreases both the crossover point v c and the 
plateau value a(v c ), but leaves the slope of the function unchanged. 



The linear dependance on v is an interesting feature. It implies that the effect of disorder 
- a measure of which is the value of o - is smaller for small v (that is, further away 
from the diagonal in the discrete form of <J a b). As we will see in the next section, the 
large-distance director correlations are controlled by the small-t> solution. Therefore 
RSB predicts that disorder has a milder effect in large-distance correlations, in marked 
contrast to the RS solution where a has a constant value. 

The full RSB form of a(v) incorporates both the linear and the constant (a' = 0) 
solutions. In fact there is a crossover between the two at a value v c , where a(v) changes 
from one regime to the other: 

J 2ixK 3 v forvG[0,i; c ] 
a{V> 1 2ixK 3 v c = const for v 6 [v c , 1] ' 1 ' 

as sketched in Fig. This form of branched function appears in spin glasses |2H] as well 
as other disorder systems where the GVM has been used, such as random manifolds [22] 
and flux lattices jHj- 

The reason that both solutions are accepted is illustrated in Fig. |21 As we shall see 
shortly, in the Eq. (J28|) . v c reduces to zero in the absence of disorder (r = 0). This in 
turn means that a(v) is identically zero for the whole range of v G [0, 1]. The no-disorder 
limit is an essential part of an acceptable solution for o~{v) and would be absent if only 
the linear branch of Eq. E3 was accepted. 

The only undetermined parameter is the crossover boundary v c , to which we turn 
our attention to. Setting v equal to v c in Eqs. (J22J) and p7|l and observing that they 
are equal, we find: 



AT 

a(v c ) = 2 7i K 3 v c = - — — exp 
v ; k B T F 



^Qrnax , 2,V C —if Qn 

H tan 



ir 2 (3K it \ir(3Kv c 
The term linear in the ultraviolet momentum cutoff (short-distance, q max = 2ir/a) is 
identical to Eq. (|TH|) encountered in the Replica symmetry limit. Since (3K ~ Q 2 /a, 
this term does not diverge. The same applies for the argument of the arctangent, which 
is also approximated by 2/(Q 2 v c ), after which the equation determining v c becomes: 



2(k B T) 2 T 

7lK 3 



exp 



2v, 



c 



H tan 



ii Q 2 ix \v c Q 



(28) 
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The approximately equal sign refers to the use of K ~ Q 2 k B T/a. 



4. Results and Discussion 

4-1. Director Correlations 

This section discusses how to obtain the domain size as a function of order parameter 
and disorder strength. To define the correlation length we examine how the director 
correlations decay. As a result of the use of a single angle to parameterize n, director 
correlations are given by |H} ITTj: 

C(r) = (n(r) ■ n(0)> = (cos [0(r) - 6(0)}} oc e - B(r)/2 , 
where B(r) = 2 / G aa [1 — cos(gr)] 

Jq 

In RSB G aa is replaced by (see |22] for detail): 



1 



Kg 2 



(29) 



<? V 9 / v^c / ^g 2 + M(v c ) 

The 1/ K q 2 term which dominates at large q originates from thermal fluctuations. The 
remaining terms, which dominate at small q (large distances), show the effect of disorder. 

In order to determine the correlations, B(r) must be calculated using the random 
field contribution of Eq. (|29)1 . The g-integration is rather involved and yields: 

B(r) = ~ I Eulerr + SI (0 - CI (0 + log f (30) 



- 2(1 - r, } smb I ^ 



sinh I — 1 — cosh ( — 



2£/ \ 2 £ 



where £ 1 = nKv c . 

Eulerr « 0.58 refers to the Euler constant, S7 and CI are shorthand notations for 
the Sinhlntegral and Coshlntergal functions, respectively. Their behavior for small and 
large arguments is given by: 

SI(x) = x + 0{x 3 ) 

CI(x) = Eulerr + logo; + 0(x 2 ) ; and CI (00) = SI(oo). 

Director correlations decay differently depending how r compares with the length 
scale £, which is identified as the correlation length or domain size: 

{exp(— r /f) , for r << £ 
(r,Y l r c ( 31 ) 

(V^J ,forr»£ 

The power-law dependance was first derived by Giamarchi and LeDoussal |Hj and is 
reproduced in our work. It shows that alignment order persists over greater length 
scales than one might first guess from the Larkin argument and the work on random 
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anisotropy magnets jSlHj. One also associates the exponential decay at short distances 
of Eq. |^ which is also found in the RS case, with the constant part of o~(v) and the 
algebraic decay with the linear part. In this context the small-w solution can be taken 
to correspond to large distance correlations. 

An exact analytic solution of v c in Eq. (|28|) is clearly impossible, however its 
dependence on the order parameter Q and disorder strength T can be deduced. It 
is then straightforward to find £ = ttKv c as a function of these two parameters, which 
is done in the next two sections. 

4-2. Domain size as a function of order parameter 

In this section we present the first theoretical model that predicts successfully the 
behavior of the domain size close to the phase transition. As discussed in the 
Introduction, the Imry-Ma argument gives £ ~ K 2 /T oc Q 2 . This expression cannot 
hold close to the nematic-to-isotropic transition since we expect that £ should diverge 
as Q — > 0. Experiments measuring £(Q) directly do not exist for nematic elastomers: 
it is practically impossible to measure simultaneously the order parameter and domain 
size of such systems. However, light scattering experiments measured £(T) the domain 
size as a function of temperature ^H] • Separately the order parameter of a monodomain 
elastomer was also measured as a function of temperature [SOI- Assuming that the 
polydomain and monodomain samples have the same form of Q(T), the data are 
combined with £(T) to make a parametric plot of £(T) versus Q(T) in Fig. El As 
expected the domain size increases at small Q. 

To determine for small order parameter we must first find the functional 

dependence of v c on Q. From their definition both v c and Q are smaller than one. 
Therefore the argument of the arctangent in Eq. (|28j) is large and tan -1 — > tt/2. Hence 

2r / 4 \ 

v ^^w^ e ^[-^ +Vc ) • (32) 

When Q <C 1 the factor of exp(— 4/ttQ 2 ) ensures that v c is also much smaller than one. 
Therefore v c <C Q~ 2 and term v c that sits on the exponent is negligible. Taking into 
account that both T and K are proportional to Q 2 for small Q, one finds: 

£ = oc Q 2 e^ (33) 

As we shall see later, for real elastomers v c ~ 10~ 4 even for Q — 1, so taking the 
arctangent equal to tt/2 is a very good approximation. 

The solid line in Fig. (J3J) shows a fit to the model expression (|33p. The agreement 
between theory and experiments is good for small Q, but breaks down for Q > 0.3. 
This is attributed to the fact that K is proportional to Q 2 only for Q <C 1, therefore 
the pre-exponential factor of Q 2 arising from K 2 /T is strictly valid only for small Q. 
A similar indirect study of the £(Q) was fitted with a power law £ oc Q~ 2 jTHj. The 
Imry-Ma correlation length was considered, forcing the authors to suggest that r oc Q 6 , 
which is hard to justify. It is not reasonable to try to differentiate the two fits, since the 
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Figure 3. Combined results of two different data sets: £(T) ^H] and Q(T) |2U] to 
form a parametric graph of The fit (solid line) comes from Eq. (|33|l and is good 

for small order parameter where the theoretical prediction holds. 



data do not come from direct measurements. However our present model is theoretically 
consistent since it assumes r oc Q 2 in accordance with the idea of crosslinks providing 
a random field that couples to the local order parameter as seen in Eq. (|Tj) . 

4-3. Domain size as a function of disorder strength 

Lattice Monte Carlo simulations have been performed to study the equilibrium ordering 
in a two-dimensional nematic system with quenched random disorder |14j . The long- 
range correlation of the director orientation was found to decay as a simple exponential. 
However, there is a flaw with the simple exponential decay, since it does not predict 
the "domain size" to diverge as T — > 0, as one must expect. The correlation length 
£ itself also decays exponentially with increasing strength of the disordering field. The 
weak disorder region was never probed because the amount of computation needed to 
produce reliable results increased sharply as T was reduced. 

In order to obtain £ as a function of T we must first find v c (T). Going back to 
Eq. ()28|) . let us rename the dimensionless prefactor as 

D = M. (34, 

71 K 6 

This is a measure of the relative disorder strength. To get an order of magnitude 
estimate for D in nematic elastomers, we consider, following the detailed experimental 
comparison in 7 ~ 0.4 x ksT, po ~ 2.5 x 10 26 m~ 3 and K ~ 1CT 12 Jm _1 , which gives 
D of the order of 10 -4 . This means that for constant order parameter Q — 1 v c is small 
and the arctangent in is equal to 71/ 2 with good precision. 

To obtain a theoretical prediction of v c (D) let us use a self-consistent method in 
Eq. (|2EJ). At Q = 1 we can write 

v c = L>exp(-4/7r + <) « 0.3De Bc . (35) 
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Figure 4. Numerical evaluation of v c (D) from Eq. (|28[) for fixed order parameter is 
plotted as solid line. The dotted line shows the analytic approximation of Eq. Ij^fijl . 
The agreement between the two curves is so good that one can hardly distinguish 
between them in the limit D — » 0. 



In the first approximation the right hand-side dependance on v c is dropped, leaving 
v c ~ 0.3 D. This value of v c is placed back in the exponent of (f33|) to give the estimate 

v c ~ 0.3De 03D . (36) 

A numerical solution of Eq. (j2HJ) can be easily found without the assumption v c -Cl and 
the subsequent approximation of the arctangent term. As seen in Fig. the agreement 
of the estimate above with the numerical solution is very good for small values of D, 
which real systems are expected to have. An important feature of this plot is that it 
starts from the origin. When there is no disorder in the system (D — 0), a(v) is zero 
for all v because v c = 0. 

The domain size has a more complex functional dependance on disorder strength 
than a simple exponential. It is inversely proportional to v c , therefore since D oc V: 

£(r)oc^e-« r , (37) 

where a ~ 0.2 (f3 2 K 3 )" 1 . The computationally obtained £(r) data points of reference 
[T4*] were fitted with the function of Eq. ()37|) . The resulting fit, Fig. Elis an improvement 
of the simple exponential function and offers, for the first time, an analytical explanation 
for the exponential decay found by Yu et al |14j . Importantly the (l/r)e _Qr function 
also predicts the divergence of the domain size at vanishing disorder, an important 
improvement of the previous work. 

Nematic elastomers can potentially provide a physical system where the spin- 
glass theory can be extensively tested experimentally. Parameter V ~ 7 2 po can be, 
in principle, controlled by changing the density or the molecular nature of network 
crosslinks, an important factor in determining many physical properties of elastomers. 

4-4- Nematic-isotropic phase transition 

The last question we address here is whether the presence of quenched disorder changes 
the character of the nematic-isotropic phase transition. In the case of XY nematics 
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Figure 5. The log-linear dependance of the domain size on disorder strength. The 
data points are taken from the simulation work ^3]. The dashed line is the fit by a 
simple exponential, provided by the authors of the original simulation work. The solid 
line shows the fit with our model expression (|57|): (9.8/r) e -°- 45r - a clearly better fit. 



whose director is confined in a plane, the transition is second order because the cubic 
term in the Landau-DeGennes expansion of the free energy density is identically zero: 
tr(Q 3 ) = 0. We calculate the energy density arising from disorder as a function of order 
parameter using the model Hamiltonian of Eq. (J0J): 



PF 



_d_ 

dm 



J \ ab q / 



d 

oc — exp 
am 



m=0 



•± J> log cr 1 ] 

q J 



(38) 

m=0 



The aim is to add F to the Landau-DeGennes expression to see if there is any significant 
change on the phase transition. In the RSB matrix algebra the trace of a logarithm is 
given by [22] : 



— tr log G 



m 



o 



1 dv, 
— log 

V 



Kg 2 



Kq 2 +[a](v) 



Substituting a(v) and integrating over q: 



1 

m 



Amq 2 dq 
(2vr) 3 



tr log G" 



X z | ^ - — 
3 2 



ttKv c tan" 



1 / irKv c 



+ l^Kql 



tan 



ttKv c 

Qmax 



Qmax 

TT 2 K 2 V 2 



(39) 



(40) 



1 + 



Qmax 



One way to simplify the above expression is to recall that v c (3K/q max 1 and therefore 
the arctangent and logarithmic terms can be expanded in a Taylor series and thus obtain 
a Q-clependent expression. However, due to the singular dependence of v c on the order 
parameter (e -1 /^ 2 ), all the terms of Eq. (J4(J|) vanish close to the critical point of this 
continuous transition. Hence the disorder energy is zero as Q — ► and the transition 
remains second order, with unchanged critical behavior. This is in contrast to the case 
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Figure 6. Plots of order parameter against reduced temperature t — T/T* that 
qualitatively show the effect of disorder. Plot (a) is obtained from the Landau- 
DeGennes theory in the absence of disorder, whereas curves (b) and (c) show the 
effect of increasing disorder. 



of 3-D-director nematics where the addition of disorder changes the phase transition 
from first to second order. The overall weak effect of disorder on the average nematic 
order in the system is sketched in Fig. |H1 Nematic order is weakened away from the 
transition, but overall there are no significant changes. 

5. Conclusion 

In this paper we have developed a quantitative description for the ordering of nematics 
with quenched disorder, following the now classical work on Replica symmetry breaking 
in spin glasses. Our model is focused on a particular case where the director is confined 
in a plane (corresponding to the 3D XY model in magnets). Such nematic configuration 
occurs in many experimental arrangements where thin films are used and is especially 
relevant for nematic elastomers. Quenched impurities (chemical crosslinks in the case 
of elastomers) act as sources of disorder by providing easy anisotropy axes. The 
competition between Frank elasticity and these random sources leads to the loss of long 
range orientational order: the mesogenic units (rod-like molecules) are assumed to be 
locally well-ordered with the magnitude of the order parameter Q uniform throughout 
the system, but the orientation of the director varies on a length scale £, called the 
domain size. We have extracted the dependence of £ on Q and T, a measure of the 
strength of disorder. In addition we have checked what effects the addition of disorder 
has on the nematic-isotropic phase transition in such XY nematics. 

The functional forms of £(Q) and £(T) adequately describe the evolution of 
the domain size, reported by experiment and simulations. In particular, £(Q) ~ 
[K 2 /Y)e 2qmaxkBT ^ n2K oc Q 2 e 4 / nQ2 provides a consistent account of the divergence of £ 
as the phase transition is approached. For Q > 0.3, where the exponential function 
varies much more smoothly, the Imry-Ma estimate (£<f ~ K 2 /V) is recovered. When 
disorder is not strong, we found the domain size to diverge as £(r) ~ (K 2 /T) e~ aV . This 
form combines the exponential decay found in some simulations with the classic Imry- 
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Ma result. It is possible that this behavior can be verified experimentally by measuring 
£ for nematic elastomers with different crosslink density. 

Quenched random-anisotropy effects were found to become irrelevant as the 
continuous phase transition of the XY nematic is approached. The parameter o that 
models the random field within the GVM is identically zero as Q approaches zero. For 
this reason, the addition of impurities only affects the state of local order, the Q(T) 
plot, well below the transition temperature where Q > 0.3. In a separate paper, we 
consider that nematics with 3D director conformation (analogous the the Heisenberg 
magnet) and find that they behave in a completely different way: the presence of disorder 
transforms the first-order transition of pure systems to a continuous one. 
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